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Abstract. An overview is presented on statistical techniques for the analytic 
evaluation of integrals for non-resonant, non-resonant depleted, non-resonant 
cut-off, non-resonant sccreened, and resonant thermonuclear reaction rates. The 
techniques are based on statistical distribution theory and the theory of Meijer's 
G-function and Fox's H-function. The implementation of Meijer's G-function in 
Mathematica constituts an additional utility for analytic manipulations and nu- 
merical computation of thermonuclear reaction rate integrals. Recent results in 
the astrophysical literature related to the use of analytic thermonuclear reaction 
rates are incorporated. 
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1 Introduction 

Understanding the methods of evaluation of thermonuclear reaction rates is one 
of the most important goals of research in the field of stellar and cosmological 
nucleosynthesis. Practically all applications of fusion plasmas are controlled in 
some way or another by the theory of thermonuclear reaction rates under spe- 
cific circumstances. After several decades of effort, a systematic and complete 
theory of thermonuclear reaction rates has been constructed (Lang, 1999). One 
of the basic ideas in this regard is that the motor of irreversibility and dissipa- 
tion is the existence of reactions between individual nuclei. The latter produce 
a randomization of the energy and velocity distributions of particles. The ef- 
fect of the reactions is balanced by the flow of the particles in a macroscopically 
inhomogeneous medium. As a result of this balance, the system reaches a quasi- 
stationary state close to equilibrium, in which steady fluxes of matter, energy, 
and momentum are present. The main ideas in the following are coming from 
statistical distribution theory and the theory of generalized special functions, 
mainly in the categories of Meijer's G-function and Fox's H-function of scalar, 
vector, and matrix arguments (Mathai, 1993; Aslam Chaudhry and Zubair, 
2002). For an overview on the appliaction and historical background of such 
functions for particle distributions, sec Hegyi (1999). An interesting application 
of them in the formation of structure in the universe, see Adler and Buchert 
(1999) and Buchert at al. (1999). A fusion of mathematical and statistical tech- 
niques enabled us to evaluate thermonuclear reaction rate integrals in explicit 
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closed forms. Some of the techniques which are used will be summarized here. 
In order to explain the ideas we will start with the evaluation of an integral over 
a real scalar variable first. Let 

r°° - / 

I(z;p,n,m)=p e- pt r np e - zt ~ dt (1) 
Jo 

for > 0, $l(z) > 0, n, m positive integers, where 3?(.) denotes the real part 
of (.). A particular case of this for np = —i>,p = 1, n = l,m = 2, 

/•CO 

I(*;l,l,2)= / e~H v e- zt dt (2) 
Jo 

is a thermonuclear function associated with equilibrium distributions in reaction 
rate theory under Maxwcll-Boltzmannian approach. As can be seen from (l)that 
the usual mathematical techniques fail to obtain a closed-form representation 
of the basic integral in (1). 

2 Statistical Techniques 

Certain special functions are related to particular probability laws governing 
products of independent exponential variables. Such laws can be related to the 
underlying physical processes. The integrand in (1) is a product of integrable 
positive functions and hence it can be made into a product of statistical densities 
by normalizing them. Consider two statistically independent real scalar random 
variables x and y with the density functions fi(x) > 0, /^(y) > for < x < 
oo,0 < y < oo and fi(x) — 0, /2(y) = elsewhere. Let u — xy, the product 
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of these random variables. Then from the transformation of the variables, u = 
xy, v = x, the density g(u) of u is given by 

»(«) = I™ -Mv)f 2 (u/v)dv. (3) 
Jo v 

The integral in (3) can be made equivalent to the integral in (1) by suitably 
selecting /i and f 2 . Let 

v = pt,u = pz m ' n , h(t) = t x - n "e-\ and f 2 (t) = e^ 4 " 7 ™, (4) 

excluding the normalizing constants. Then 

f°° 1 Z" 00 1 / 

/ -fi(v)f 2 (u/v)dv = / -e-^{ V t) X - np e- zt ~ dt. 
Jo v Jo t 

Hence 



roc i roo 

p pp -e- pt {pt) 1 - np e-* t ' dt = p e pt t~npe- zt ~ dt 
Jo t Jo 

which is exactly the integral to be evaluated in (1). We have identified the 
integral as the exact density of the product of two real scalar random variables, 
u = xy. Since this density is unique the idea is to evaluate this density through 
some other means. Notice that u is a product of positive variables and hence the 
method of moments can be exploited profitably. Consider the (s — l)th moment 
of u, denoted by expected value of u s_1 . That is, 

Eiu 3 - 1 ) = Eix^Eiy 8 " 1 ) 

due to statistical independence of x and y. Let 

gi (s) = Eix*- 1 ) and g 2 (s) = E^- 1 ). 



Then gi(s) and (72 (s) are the Mellin transforms of /1 and fi respectively. Then 
from the inverse Mellin transform the unique density of u is available as 

g(u) = J u- s g 1 (s)g 2 (s)ds,i = \/ z T. (5) 

where L is a suitable contour. But, excluding the normalizing constants, 
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(s) = / i-ViWd* 
Jo 

= / fl 



and 



= r(l - np + s) for 5R(1 - np + s) > 0, . (6) 

52(5) = / t s -\f 2 (t)dt 

Jo 

oc 





777 

— r(ma/n) for > 0. (7) 



Then 



77^- / «- s .gi( S )g 2 (s)d S - /(m/n)r(m S /n)r(l-n /0 + .s)(pz m /")- s d S . (8) 

27TI 7 L 27TI 7 L 



Therefore 



P / °° e-PH- n Pe- zt ~ n/m dt 



y -e- p *(pf) 1 -"''e-^" r ' /m di 

p«P( TO / n )_L /* r(ms/n)r(l -n/7 + s)(pz m /™)- s d S 
2tt« 7 l 

mp np — ( r(ms)T(l-np + ns)(z m p n )- s ds (9) 
2tt« J Li 



by replacing s/n by s. Our aim now is to evaluate the contour integral on the 
right side of (9) explicitly into computable forms. The contour integral in (9) 
can be written as an H-function which can then be reduced to a G-function 
since m and n are positive integers. 

3 G- and H-Functions 

For the sake of completeness we will define G- and H-functions here. For the 
theory and applications of G-functions see Mathai (1993) and their relation to 
incomplete gamma functions (Aslam Chaudhry and Zubair, 2002). Let <p(s) and 
ip(s) be the following gamma products. 

{n^i£fe + M) {n^rq - aj - a jS )} 

ns) {nj =m+1 r(i - b 3 fa)} {TP j=n+1 T( aj + a jS )} 

and 

(f)(s) = il>(s) with a 3 = l,j = l,...,p,p k = l,k=l,...,q. 
Then Fox's H-function, denoted by 

rrm,nj \ _ frm,nr 

n p,q \ Z ) ~ n p,q L Z I(6 1 ,/3 1 ),...,(6 9 A) J 

is defined as 

H£f(z) = f 4,{ S )z- s ds (10) 

JL-2 

where L 2 is a suitable contour separating the poles of the gammas T(l — aj — 
(Xjs),j = 1, ... ,n from those of T(bj + (3js),j = 1, . . . , m. The theory and 
applications of H-functions are available from Mathai and Saxena (1978). In 
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(10) when aj = l,j = l,...,p and ft = l,fc = we have Meijer's 

G-function denoted by 



U p,9 — ^p,q i Z \b 1 ,...,b q J 



= f <j)(s)z- s ds (11) 



>L 3 

where L 3 is a suitable contour. Various types of contours, existence conditions 
and properties of the G-function are available from Mathai (1993). 

Now, comparing (9) and (10) our starting integral is evaluated as follows 

I(z;p,n,m) = mp np - — / T(ms)T(l — np + ns)(z m p n )~ s ds 

27TZ J Ll 

= mp n "H^[( Z m p n )\^ mU1 _ np!n) ]. (12) 

The H-function in (12) can be reduced to a G-function which can again be 
reduced to computable series forms. For this purpose we expand the gammas in 
the integrand in (12) by using the multiplication formula for gamma functions, 
namely, 

T(mz) = (2n)^m mz 'ir(z)T{z + . . .T(z + ^— -) (13) 
m = 1,2, . . . 

Expanding T(ms)T(l — np + ns) by using (13) we have 

mp np T(ms)T{l - np + ns) = {2n)^ 2 - m -^p np min- np+ i (m m n n ) s 

1 711 — 1 

x T( S )T(s + -)...T(s+ ) 

m m 

1 — np s , n — np. 

x T(s+ -)...T(s + f-) 

n n 



Substituting these back and writing as a G-function we have 

I{z;p,n,m) = p np {2^ {2 - m - n) m^ n - np+ ^ 

m n 

x ^[^lo,*,..,^,^,..,^]- (14) 

4 Non-Resonant Thermonuclear Reaction Rate 

The Maxwell-Boltzmannian form of the collision probability integral for non- 
resonant thermonuclear reactions is (Lang 1999; Haubold and Mathai, 1984; 
Hussein and Pato 1997; Bergstroem et al. 1999; Ueda et al. 2000) 

r oo 

h = 7(z;l,l,2) = / y v e-ye- z y- 1,2 dy 
Jo 

= ^G^\ ^ 1+v \. (15) 

In stellar fusion plasmas the energies of the moving nuclei are assumed to 
be described by a Maxwell-Boltzmann distribution, Eexp(—E/kt), where T is 
the local temperature and k the Boltzmann constant. Folding the cross section 
of a nuclear reaction, a(E), with this energy (or velocity) distribution leads to 
the nuclear reaction rate per pair of nuclei: 

poo 

< av >= {8/nn) 1/2 (kT)- 3/2 / <r(E)exp(-E/kT)dE, 

Jo 

where v is the relative velocity of the pair of nuclei, E is the center-of-mass 
energy, and /j, = TO1TO2/ (mi + m-i) is the reduced mass of the entrance channel 
of the reaction. In order to cover the different evolution phases of the stars, i.e. 
from main sequence stars to supernovae, one must know the reaction rates over 
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a wide range of temperatures, which in turn requires the availability of a(E) 
data over a wide range of energies. I is the challenge to the experimentalist to 
make precise <r(E) measurements over a wide range of energies. For the class of 
charged-particle-induced reactions, there is a repulsive Coulomb barrier in the 
entrance channel of height E c = ZiZ 2 e 2 /r, where Z\ and Z 2 are the integral 
nuclear charges of the interacting particles, e is the unit of electric charge, and 
r is the nuclear interaction radius. Due to the tunneling effect through the 
Coulomb barrier, a(E) drops nearly exponentially with decreasing energy: 

a(E) = S(E)E- 1 exp(-27rr 1 ), 

where r\ = ZxZ^e 1 jhv is the Sommerfeld parameter, h is the Planck constant. 
The function S(E) contains all the strictly nuclear effects, and is referred to as 
the astrophysical S(E) factor. If the above equation for a(E) is inserted in the 
above equation for the nuclear reaction rate < av >, one obtains 

/>oo 

< av >= (8/^) 1/2 (fc7T 3 / 2 / S(E)exp(-E/kT - b/E 1 ' 2 )dE, 

with b = 2(2^) 1 / 2 7r 2 e 2 ZiZ 2 //i. Since for nonresonant reactions S(E) varies 
slowly with energy, the steep energy dependence of the integrand in the equa- 
tion for < av > is governed by the exponential term, which is characterised 
by the peak near an energy EO that is usually much larger than kT, the mean 
termal energy of the fusion plasma. The peak is referred to as the Gamow 
peak. For a constant S(E) value over the energy region of the peak, one finds 
E = (bkT/2) 2 / 3 . This is the effective mean energy for a givern reaction at a 
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given temperature. If one approximates the peak by a Gaussian function, one 
finds an effective width S = A(EokT) 1 / 2 /3 1 / 2 . In the following, this approxi- 
mation is not made and the respective integrals are analytically represented, 
beginning with 

<av> = ( — )2 > —j- ^ 

x / e- y y v e- zy dy 
Jo 

where S^fi) denotes the /z-the derivative. The G-function in (15) can be ex- 
pressed as a computable power series as well as in closed-forms by using residue 
calculus. Writing the G-function in (15) as a Mellin-Barnes integral we have 

^[^| ,i,i+,] = ^- J r»r(.s + i/2)r(i + „ + s )(z 2 /4)-^ s . (ie) 

Case (1): v ^ A = 0, 1, 2, . . . In this case the poles of the integrand are 
simple and the poles are at the points 



a = 0,-1,-2,...; 




= -1 - z/,-2- v, . .. 



1 7 2 

E ^ r(--r)r(i + i/-r)( T r 

r=0 

= r(i)r(i + I /) F 2 (;i -i/;-^) (17) 

where, in general, p F q denotes a general hypergeometric function. The sum of 
the residues at s = — |, — | — 1, ... is 

oo i 1 2 

E ^ r(-i-r)r(i + „-r)(fL)**- 

r=0 
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n-\)n\ + v)C^) h »f^\,\-v- z -). (is) 



The sum of the residues at s = — 1 — v, — 1 — v — 1 , . . . is 



00 1 2 

£ ^ r(-i-,-r)r(-i-,-r)(^)^ 

r=0 



= T{-l-v)Y{- 1 --u){ Z -)^ F 2 {;y + 2,u+\-- Z -) (19) 

Then from (15) to (19) we have 

7(z;l,l,2) = / y v e- y e- zyl/2 dy 
Jo 

= """^O^I^jHo.i.l+L/] 

= r(l + z/) F 2 (;i 

- 2r ( - + ,)( T ), F 2( ;-;--,;- T ) 

r(-i-^)r(-i-^) ,2 

7(1) ( 4 j 

3 z 2 

X f 2 (; V + 2,1/ +-;- — ) (20) 

for v 7^ A = 0, 1, 2, . . . When v is a positive integer the poles at s = — 1 — 
nu, — 1 — i/ — 1, . . . are of order 2 each. Hence the corresponding sum of residues 
can be written in terms of psi functions. Similarly when v is a negative integer, 
positive or negative half integer the corresponding sums will contain psi func- 
tions. Details of the computations can be seen from Haubold and Mathai (1984) 
and their application by Hussein and Pato (1997), Bergstroem et al. (1999), and 
Ueda ct al. (2000). 
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5 Modified Non-Resonant Thermonuclear Re- 
action Rate: Depletion 

With deviations from the Maxwell-Boltzmann velocity distribution of nuclei in 
the fusion plasma, a modification which results in the depletion of the tail is 
introduced in Haubold and Mathai (1986a); see also Kaniadakis et al. (1997), 
Kaniadakis ct al. (1998), and Coraddu et al. (1999). In this case the collision 
probability integral will be of the following form: 

poo 

h= y»e- v e- zv ~ 1/2 dy. (21) 
Jo 

We consider a general integral of this type. Let 

I(z;S,a,b,m,n) = ^e" *" 6 * ~ zt n dt. (22) 
Jo 

Expanding e~ htS and then with the help of (14) one can represent (22) in terms 
of a G-function as follows: 

I(z;S,a,b,m,n) = ^ S_L a -(p+kS+i) ^^^-rn-^^^+p+kS 

k=0 

X g "m+J mron „ lo,^,...,2^1,l±£±M,...,g±£hM] (23) 

for m.[z) > 0,3?(a) > 0, 5R(6) > 0, m, n = 1, 2, . . . . the case in (21) is for a = 
1, b = 1, n = 1, m = 2. With v = p + kS the G-function in (21) corresponds to 
that in (15). When S is irrational and p is rational, the poles of the integrand will 
be simple and the G-function is available in terms of hypergeometric functions. 
Other situations will involve psi functions. From the asymptotic behavior of the 

12 



G-function, see for example Mathai (1993), one can write the integral in (21), 
for large values of z as follows: 

h « ^{P/i) 2 -^e-^l^\p = 3(z/2) 2 / 3 . (24) 

6 Modified Non-Resonant Thermonuclear Re- 
action Rate: Cut-Off 

Another modification can be made by acut-off of the high-energy tail of the 
Maxwell-Boltzmann distribution. In this case the collision probability integral 
to be evaluated is of the form 

r d 

h= t- p e- at - zt dt,d <oo, (25) 
Jo 

see the details from Haubold and Mathai (1986b). We will consider a general 
integral of the form 

f d 

I(z;d,a,p,n,m)= y- np e- av ~ zv ~" /m dy, d < oo. (26) 
Jo 

In order to evaluate (26) explicitly we will use statistical techniques as discussed 
earlier. Let x and y be two statistically independent real random variables 
having the densities Ci/i(x), < x < d and 02/2(2/), < V < 00 with fi{x) and 
f2{y) equal to zero elsewhere, where c\ and C2 are normalizing constants. Then 
taking 

h{x) = x- n P +1 e- ax and f 2 (y) = 
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and proceeding as before one has the following result: 
I(z;d,a,p,n,m) = / / "' ' ; <U 



m5 n -i(2^)( 1 -™)/ 2 rf-»P+ 1 ^ 



(-ad) r 



r=0 

T+2+j-l 



d n m m l-p+ r+1 +^- 1 ,j=l,...,n;^,j=l,. 



(27) 



for 3?(z) > 0, > 0, SR(o) > 0. Then 

h = I t- p e- at - zt dt 
Jo 



x G'l 



3,0 



fi-jO+r+2 
4rfl-p+r+l,0,| 



(28) 



For large values of z the G- function behaves like n^x 2 e 2x2 ,x = |g , see for 
example Mathai (1993). Then for large values of z, 



T * ™ d ~ p+1 (j^) 2 e -° d - 2(22/4d)S . 



(29) 



Explicit series forms can be obtained for various values of the parameters with 
the help of residue calculus. For example, for $t(z) > 0, d > 0, 9?(a) > 0, — p + 
r + 1 ^ fj,,/j, = 0,1,... 

f d 

h = I t- p e~ at - zt dt 
(-ad)' 



7r-id-o +1 Y2 



~ (-i)T(i-^) z \ 

^ i/!(-p + r+l-!/) l 4rf j 



14 



where 



A = 4>{(i + 1) + 4>{-n + -), 4>(z) = ^lnT(z) 



l 



Other cases and further details may be seen from Haubold and Mathai (1986b). 

7 Screened Non-Resonant Thermonuclear Re- 
action Rate 

Taking into account plasma corrections to fusion processes due to a static or 
dynamic potential, see the details from Haubold and Mathai (1986c), Shaviv 
and Shaviv (2001), Wierling et al. (2001), Bahcall ct al. (2002), the collision 
probability integral to be evaluated is of the form 




(31) 



We will consider a more general integral of the following form: 




oo 



I(z; t; a, v, n, m 



,v e -ay-z(y+t) n 



m 



dy. 



(32) 



This can be evaluated with the help of (14) and (27). Let 



Ni(z; a, v, n, m) 




15 



Then from (14) 



— m — l n + u n + u 



m"'n" 



(33) 



for Sft(a) > 0,3t(z) > 0,m,n positive integers. Let 

N 2 (z) — N 2 (z;d,a,v,n,m 
,-d 



- f 

Jo 



Then from (27) 



i^+i {-adf 



N 2 {z) = {2^ 1 ~ m ^ 2 m^n- 1 d v+1 Y^ 



x G 



r! 

, + r+2+j-l j = 1 )W 



m+n,0 r ^ 

n,m+n[^ nmm I ,■,•=!,...,„, ^ ,j=l,...,m J 



(34) 



for 9?(z) > 0, d > 0, 3?(a) > 0, m, n positive integers. Then by a change of 
variables and rewriting the integral we have 



v -ay-z(y+t) n/r ' 



V e 



dy 



f +1 e 0l £(r)(-l) r 



r=0 



x [JVi (zi; ai, ^ — r,n, m) — N 2 {z\; 1, ai, ^ — r, n,m)] (35) 



where 



ai = at, zi = zt, ) 



,0! = 1, 



r\{v — r)\ ' 

and iVi(.) and N 2 {.) are defined in (33) and (34), respectively. 
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8 Resonant Thermonuclear Reaction Rate 

When a resonance occurs in the low-energy range of nuclear reactions the col- 
lision probability integral can have the following form, for details see Haubold 
and Mathai (1986d). 

R(q,a,b,9) = j o t" (& _ y? - ~ 2 dy. (36) 
We will consider a more general integral of the following type: 

/-co e -at-qt- n/m 

R 1 (q,a,b,g,v,n,m)= t" '— — — -dt. (37) 
Jo (o - *) + 9 

The method to be employed here for evaluating this integral is to replace the 
denominator by an equivalent integral. That is, for g 2 > 0, 

1 



(b-t) 2 +g 2 

Then expand 



roo 

/ e-W-'f+ft'dx. (38) 
Jo 



k=0 ' fci=0 

Substituting back in (37) , integrating out t first and then integrating over x one 
has the following result: 

/>OC 

/ t v+kl e- at - qt ~" /m dt 
Jo 



^m+n,0 
A ^0,m+n 



q m a n 



^loi,...,H^i,i±^l,...,^±i 



(39) 



and 



f°° 2 A-' 
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From (39) and (40) 



Ri (<?, a,b,g,v,n,m) 
(-1) 



k ->' 
fi 

k=0 w ' fci=0 



2fc-fci 



x a 



(,y+l + fei) ( -2 7r )i(2-n-m) TT1 5 T1 ^+ fe i+5 



^m+n,0 



m" l n 



m 2 n" 



-1 l+f+fci n+K+fej 



Then for n = 1, m = 2, i/ = 0, 
i?(g,a,6,g) 



i 2„ ^ 



(~l) fc 
-'5 2 « (5 2 ) fc 

x E (2?) ^jH 2 * - * 1 

fei=0 



x t7 3 ' 

x u 0,3 



r 2 

fa 



4 IU <2 



o,i, fci 



for (5z£L<i )I/= 



(41) 



9 Computations 

For computational purposes we will consider the four basic integrals associated 
with the cases: non-resonant reactions, non-resonant "cut-off" reactions, non- 
resonant screened reactions, and non-resonant "depleted" reactions. Let 

Mz,u) = / y»e-y- z y- 1/2 dy 
Jo 

Mz,d,v) = [ y"e-y- z y- 1/2 dy 
Jo 

r°° _ i 

Mz,t,u) = / jfe-y-'b+Q ^dy 
Jo 

poo 

Mz,S,b,u) - / y v e- v - b «°- xv ~ 1 "dy. (42) 
Jo 
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The exact expressions for these are given in (15), (28), (35), (21) respectively. 
The symbolic evaluation of all these integrals can not yet be achieved with Math- 
ematica (Wolfram, 1999). Those inetrals that involve no singularities are done 
by taking limits of the indefinite integrals. The definite versions of the integrals 
are done using the Marichev-Adamchik Mellin transform methods (Adamchik 
1996). The integration results are initially expressed in terms of Meijer's G- 
function, which are subsequently converted into hypergeometric functions using 
Slater's theorem. The notation for Meijer's G-function, belonging to the imple- 
mented special functions of Mathematica, is 

MeijerG[{{ai, a„} , {a n+1 , a p }} , {{&i, b m } , {b m+1 , ...,b q }} , z]. (43) 

Analytic expressions for the following Meijer's G- functions are available on Wol- 
fram Research's Mathematical Functions web page (Wolfram, 2002): 

G{m,n,p, q} = G {3, 0,0, 3} = http://functions.wolfram.com/07.34.03.0948.01, 

(44) 

and 

G{m,n,p, q} = G {3, 0,1,3} = http://functions.wolfram.com/07.34.03.0955.01. 

(45) 

Approximations for large values of z can be worked out with the help of the 
asymptotic behavior of G-functions, see for example Mathai (1993). These are 
the following for z very large: 
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3( Z 2 /4)*-6( 2 2 /4) 5/3 



(46) 



A number of exact and approximate graphs of these four integrals, for vari- 
ous values of the parameters, using the integration routines in Mathematica 
(Wolfram, 1999), are available from Anderson, Haubold and Mathai (1994). 

10 A Generalization 

A mathematically interesting integral corresponding to (1) can be evaluated. 
Consider the integral 



Then take fi{x) = axPe- px ,x > 0, f 2 (y) = c 2 e- y \j > 0,y > and fx(x) = 
0, f 2 {y) = elsewhere, where c\ and c 2 are normalizing constants. Then u = 
xy = z 1 / 7 and from (3) one has the integral in (44) evaluated as the following: 

/ = (w P r 1 H 2 ^[pz 1 ^\ (pAh{aAh) },0< z < co, 

where H(.) is the H-function defined in (10). When 7 is rational the H- function 
can be rewritten in terms of a Meijer's G- function and then (44) can be evaluated 
in terms of the result given in (14). For specified values of 7 and p one can obtain 
computable representations for the H-function. 




(47) 
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